namespace Likelihoods{

class Logit : public Density::GeomBridge
{
		friend class Density::GeomBridge;
	public:
		Logit(double *Y, mat X, Distribution::Distribution *S,mat s) : Density::GeomBridge(Y, X, S){ 
			_s=s;
			_is=inv(s);
			//d= new boost::math::normal_distribution<>(0,1); 
  }
		~Logit()
		{
			//delete d;
		}
		inline double Likelihood(mat theta1)
		{
			mat theta=theta1.t();
			double L=0;
			double sum=0;
			int m=theta.n_cols;
			//Xb
			int n=Get_n();
			mat X=Get_X();
			double *Y=Get_Y();
			for(int i=0;i<n;i++)
			{	
			//	cout << theta.col(0) << "\\";
				mat foo= X.row(i);
				sum=dot(foo,theta);
				double phi=Phi(sum);
				double phi2=1-Phi(sum);
				if(Y[i]==1){
					L+=log(phi);
				}else if(Y[i]==0){
					L+=log(phi2);
				}else{
					cout << "Warning values of Y!= {0,1}";
				}
			}
			//cout << exp(L);
			return L;

		}
		inline double Phi(double x)
		{
			double res=1/(1+exp(-x));
			if(res==0)
			{
				res=0.00000001;
			}else if(res==1){
				res=0.99999999;
			}
			return res;

			//return boost::math::cdf<double, boost::math::policies::policy<> >(boost::math::normal_distribution(0,1),x);
		}

		double Prior(mat theta){
			double res=0;
			mat foo=theta.t();
			int d=theta.n_elem;
			mat bar=foo*_is*foo.t();
			//cout << "//" << det(_s);
			double de=det(_s);
			if(de>100000){de=100000;}
			res=-d*0.5*(log(2*PI))-0.5*log(de)-0.5*as_scalar(bar);
			return res;
		}
		mat GradLik(mat theta)
		{
			mat X=Get_X();
			double *Y=Get_Y();
			int p=X.n_cols;
			mat sum(p,1);
			sum.fill(0);
			int n=X.n_rows;
			double xb=0;
			for(int i=0;i<n;i++)
			{
				xb=dot(X(i,span::all),theta);
				if(xb!=xb)
				{
					cout << "cb3";
				}
				double gpxb=Phi(xb);
                                double pxb=exp(-xb);
				
				//cout << "foo1 "<<foo1;
                                //cout << "foo2 "<<foo2;
				mat temp=gpxb*pxb*X(i,span::all);
				if(Y[i]==1)
				{	
					sum=sum+temp.t();				
				}else{
								
					sum=sum+temp.t()+X(i,span::all).t();				
				}
			}
			sum+=-as_scalar(arma::sum(theta,0))/_s(0,0);
			for(int i=0;i<p;i++)
			{
				if(sum(i,0)<-1000){
					sum(i,0)=-1000;
				}else if(sum(i,0)>1000){
					sum(i,0)=1000;
				}
				if(sum(i,0)!=sum(i,0))
				{
					cout << "2";
					cout << dot(X(i,span::all),theta) << "\n";
				}

			}
			
			return sum;

		}
	        mat Get_s(void){return _s;}	
	private:
		mat _s;
		mat _is;
		//boost::math::normal_distribution<> *d;
	/*private:
		double _phi_n;
		double _phi_n1;
		mat _X;
		double *_Y;
		int _p;
		int _b;
		int _n;*/
};

} // end namespace

